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ABSTRACT 



This study is to understand the effects of nonlinear 
behavior of structures using a finite element method. 

The nonlinear behavior equations are derived from 
equations of motion and constitutive equations. The basic 
theory is the principle of virtual work. 

The thesis begins with a comprehensive formulation of 
continuum-based finite element theory. The theoretical 
portion concludes with the details of both spatial and 
temporal discretization, including a discussion of 
nonlinearity . 

In particular, large displacement problems and visco- 
elastic problems remain a challenging engineering problems 
today. The viscoelastic problems depend on the relaxation 
function which is the source of material nonlinearities. 
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I. INTRODUCTION 



A. ANALYSIS OF STRUCTURES 

Analysis of structures subjected to large displacements 
requires the inclusion of nonlinear terms in the strain 
displacement relations. In the formulation of the problem, 
either total Lagrangian or Eulerian approach is adopted. In 
order to assure structural integrity for loads beyond linear 
range, nonlinear analysis is very much mandated. 

Another class of problems is the viscoelastic behavior 
of structures - problems where stresses in the components 
vary with time at a constant deformation state or strains 
vary with time at constant stress states. 

These problems as applied to simple one-dimensional 
structures will be addressed in this research work. In what 
follows is a brief survey of literature in this area. At 
the end of this chapter, overview of the thesis is 
described. 

B. LITERATURE REVIEW 

A method for nonlinear static and dynamic analysis using 
cable and beam elements is given by Bergan et al [1]. 
Formulation allows for material nonlinearities and very 
large deformations. Gadala and Oravas [2] present a 
comparison of finite element formulations based on the 
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Lagrangian, updated Lagrangian and Eulerian formulations for 
nonlinear mechanics problems. 

Henriksen [3] applies a single integral constitutive law 
to describe viscoelastic behavior in his formulation. 
Ghoneim et al [4] treat the total strain rate as sum of 
viscoelastic and viscoplastic contributions in implementing 
the constitutive equations in their formulation, Pinsky and 
Kim [5] use variational approach for studying the nonlinear 
viscoelastic shell behavior by representing the layered 
shell as a multi -director field. The formulation includes, 
within individual layers, the effects of transverse shear 
and transverse normal strain to arbitrary orders in the 
layer thickness coordinate [6, 7], 

Saigel and Yang present a quadrilateral shell element 
and describe implementation of elastic-viscoplatic 
constitutive law for large deformation. A large collection 
of papers on the application of ADINA may be found in [8] . 

Response of nonlinear mechanical systems to random 
excitation is given by To [9] . Thermomechanical response of 
uniaxial bars with thermoviscoplastic constitutive laws is 
discussed by Allen [10]. Large deformations of the 
contacting bodies on the contact surfaces and solution 
methods are given consideration in [11, 12]. Large 
deformations of composite beams subjected to axial dynamic 
load is reported in [13]. 

Based on Eulerian-Lagrangian formulation, frictional 
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contact analysis is discussed by Haber et al [14]. 

Halleux and Casadei [6] present transient large strain 
analysis of structures based on biquadratic elements. 
Polar-decomposition is used in arriving at explicitly 
derived tangent stiffness matrices for a frame element by 
Kondoh and Atluri [15]. A method based on co-rotational 
coordinate system, where the local cartesian coordinate 
system attached to each element translate and rotate is 
described by Mattiasson [16]. 

Inelastic problems based on incremental virtual work 
equations is given by Nagtegaal [17]. A constant strain 
rate is assumed during an increment in integrating the 
constitutive equations . 

Planar curved beams undergoing large rotations is 
analyzed using Helinger-Reissner mixed variational principle 
by Noor [18]. Steele discusses several nonlinearities 
including buckling, creep, and nonlinear elastic problems in 
[19]. The use of Lagrangian formulation using an elasto- 
plastic model is presented by Voyiadjis and Navaee [20]. 
Hirashima et al [21] present an approach for the snap- 
through analysis of curved shell panels and determination of 
initial stress matrix. 

Importance of viscoelastic relaxation effects on a 
pellet are investigated by Basombrio [22] , where it is 
observed that the duration of the ramp input influences the 
response considerably. Viscoplastic constitutive models are 
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obtained explicitly, without resorting to inversion, to be 
directly implemented in a finite element formulation [23] . 
An intrinsic time scale that can vary with time and position 



is used, together with means to control the numerical 
integration in the approximation of constitutive model by 
Chambers and Becker [24], 

Chandra et al [25, 26] discuss applications of visco- 
plastic models in the metal -forming . Elasto-plastic models 
are used in analyzing unilateral contact and friction 
problems by Cheng et al [12]. Chung et al [27] suggest 
methods for varying parameters in guiding solutions for the 
simulation of plastic forming processes. Liu et al [28] 
propose formulation for path-dependent materials based on 
arbitrary Lagrangian Euler method and include an elastic- 
plastic wave propagation problem. 

Rheinboldt and Riks describe several algorithms based on 
continuation method and mesh refinements for the solution of 
nonlinear finite element equations. Ryu and Arora [29] 
propose an approach to solve nonlinear finite element 
equations based on substructures concept. 

Several papers pertaining to visco-elastic, visco- 
plastic behavior and their implementation is reported in 
[10, 7, 4, 30, 31]. 
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C. THESIS OVERVIEW 



The objective of this thesis is to review some of the 
literature available in the analysis of large deformation of 
structures, to develop the theory to include large 
displacements and viscoelastic behavior of one dimensional 
structures, and in the process, to understand the physics of 
the effects of nonlinearities. Chapter II presents 
theoretical formulation based on principle of virtual work 
and incremental load method. Descretizat ion based on a 
simplex element is given. 

Chapter III gives the algorithms used in the solution 
procedure. Full Newton-Raphson and Modified Newton-Raphson 
methods are described. 

Chapter IV describes some test cases and numerical 
results obtained. 

Finally, Chapter V provides conclusions and some remarks 
for future work in this area. 
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II. THEORETICAL FORMULATION 



A. OVERVIEW 

In this chapter we consider a nonlinear continuum 
equations for a rod element. This approach is based on the 
principle of virtual work [48]. The resulting equations are 
in convenient form for numerical solution via finite element 
methods [49]. This, of course, is the conceptual basis of 
the discretized-solid finite elements, introduced by 
J.T.Oden et al [50] in the elastic structures and nonlinear 
continua . 

B. NONLINEAR CONTINUUM EQUATIONS 

Discuss a rod element with the total Lagrangian Method. 
It may be proved that internal virtual work is equal to 
external virtual work. Namely a deformable system is in 
equilibrium if and only if the total external virtual work 
is equal to the total internal virtual work for every 

IS' 

virtual displacement consistent with the constraints. 

1 . Principle of Virtual Work 

In this section we prove that total internal virtual 
work is equal to total external virtual work. Consider the 
concept of minimum potential energy [51]. The condition 
is satisfied to equilibrium state and the change of 

potential with respect to the displacement 
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must remain stationary. 



And 



consider the load- 



displacement and stress-strain relationship of figures 

shown below. 




-^sul 4 - 



Figure 2.1 Load-Displacement Curve 




Figure 2.2 Stress-Strain Curve 



In General, the total potential energy is equal to the 

summation of strain energy and potential energy. 

ie , 7tp = U + W (2.1) 

where, 7tp : Total potential energy 

U : Strain energy (internal work) 

W : Potential energy of load system (external 
work ) 
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From Figure (2.1), 

W = - J p 6u 
= - 2 Pi Ui 

where, Pi : it h load 

ui : deflection under i th load 

From Figure (2.2), 

6 

U = J I a SE d(vol) 

VO 1 o 



( 2 . 2 ) 



= .5 



a € d(vol) 



V O 1 



(2.3) 



Substituting the equation (2.2) and (2.3) in equation (2.1), 
we get 



. 5 a € d( vol ) - 2 p u 

vo 1 i i 



(2.4) 



Equation (2.4) can be represented by matrix form, the 
following is the matrix formula. 



7Up 



where , 



! .5 {a} {€} d(vol) - { u. > { p. } 



(2.5) 



{a} 
{€} 
{Ui } 
{Pi > 



a vector of stresses 
a corresponding vector of strains 
a vector of elemental nodal displacements 
a vector of elemental nodel forces 



Let us define [B] , a matrix relating strains and 
displacements , 



ie, {€} = [B] {u> (2.6) 

and define [D] , a matrix relating stresses and strains, 

ie, {a} = [D] {€} (2.7) 
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Substituting equation (2.6) & (2.7) in equation (2.5), we 

get 






T T 



7Tp = 

Now we 
is as 



J (ui } 

V O 1 

consider 
follows : 



{B} [D] [B] d(vol) {ui } 

the minimum potential 



T 

- {Ui } 

energy, 



(Pi } (2.8) 
the condition 



d?rp 

o 

d{ui } 

Taking derivatives with respect to (ui } 
and applying to equation (2.9), we get 

{pi } = J [B] CD] [B] d(vol) {ui } 

V O 1 



(2.9) 

in equation (2.8) 



( 2 . 10 ) 



From relating loads and displacements, we can get the 
elemental stiffness matrix [k] . 



ie, [K] = J [B] [D] [B] d(vol) (2.11) 

V O 1 

Summing all the elements of the entire structure, equation 
( 2 . 10 ) becomes 

{P} = [K] {u} (2.12) 

where, {P} =2 {Pi > 

[K] = 2 [K] 

[u] = 2 {ui> 

2 . One Dimensional Finite Element 

Consider a shape function which is linearly function 
with geometric boundary condition with one dimensional 
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finite element. We define that the shape function, 



'Pn , is 



( x ) = a N + b N X 



In order to derive the constants in the shape function for a 
rod element with linear displacement field, use the 
geometric boundary condition. 

We may obtain the constants as follows: 

X2 

a = 

i Lo 

1 

b = (2.13) 

i Lo 



2 Lo 

1 

b = 

2 Lo 



The relationship between displacement and shape function is 
given below 
M 



u = n !^ n u n (N=l; M=2, . . . . ) 

= 'J 7 ! ( X ) Uj ( t ) + ^ ( x ) u 2 ( t ) 



(2.14) 



Another consideration is the viscoelastic materials. For 
viscoelastic materials, the stresses may be expressed in 
terms of relaxation functions as given by 



S ( t ) - i 



G(t - £ ) s e (f ) d£ 
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I(t) 



( 2 . 15 ) 



= S (t) - 



where , 



e e 

S ( t) = S [r (t ) ] 



I(t) 



u 

O 



G ( t - | ) 

~dT~ 



: elastic response 

Se ( £ ) df 



G(t): relaxation function. 

Assume the element is in equilibrium and quasistatic state. 
The motion equation for a one dimensional finite element is 
given by J . T . ODEN [ref. 6], 

J S\^ ( 1 u ) dV = P (u ,t) 

Vo N, x M, x M o N M 



where, Vo = Ao Lo : undeformed volume 

P : total generalized force at node N 

N 



On substituting for above shape functions and their 
derivatives, we may obtain 

f N (u M , t, I) = Vob N [ S e - I] (1 + b k u k ) - P N (u M ,t) 

=0 ( 2 . 16 ) 



C. CONSISTENT NODAL FORCES 

Consider a rod element in two configurations: undeformed 
and deformed. 

Let the distributed forces per unit length in these 
configurations be q and q* respectively. The consistent 
load vector depends on both the magnitude of the load and 
the deformed state. 
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One dimensional rod element under distributed load is 



showed as following figures. 



. „ __ %* c * , 



i a - 




• \ 




2, ■■ 


u 


1 1 


k- - 




4 



U 



where, Lot undeformed rod element length 
L : deformed rod element length 

q(x,t) is distributed load per unit undeformed length 
q*(y,t) is distributed load per unit deformed length 

Figure 2.3 One-Dimensional Rod Element under Distributed 
Load . 



Consider the relationship between the deformed and 
undeformed coordinate. It may be written as 

y(x,t)-x+u(x,t) (2.17) 

where x: undeformed coordinate of a point 
y: deformed coordinate of a point 
u: displacement of a point within element 
The consistent nodal forces vector is given by 

P = l PFN dV -- 1 PF 0 NdV. 

In index notation, nodal forces may be written as 

P = fp^dV = f P.F,* dV„ 

N Jv N i8 7v 0 N 



(2.19) 



with 



N = 1, 2 

where F : body forces per unit mass deformed 

Fo : body forces per unit mass undeformed 
P : deformed density 

P undeformed density 

o 

Relating q*(y,t) and q(x,t) is given by 
q( x , t ) = \ q* ( y , t ) 
where, \ = dy/dx : stretch 

Consider the undeformed body force per unit mass, 
may write 

p- . Vda _ ^ x 

r «- dm - ^ 0 a 0 ' p. A „ 

Using equation (2.17) 

y = x + u 

M M 



dy 

= 1 + \|/' u 

dx M, x M 



The undeformed body force per unit mass is then 

q* 

o o Ao 



A 

F = 



( 1 +\|/' M, x UM ) 



The consistent nodal forces take the following form, 

p n = L q* + I 



M^N dx 



Substituting the shape function in equation (2.25), 
get 



(2.20) 

then we 
( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 

(2.24) 

(2.25) 
we may 
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q*(a N + b N x) dx 



(2.26) 



= (1 + 



b M U M 



I 



Lo 



Now we consider some special cases of equation (2.26). 
One case is the uniformly distributed loads. 




Figure 2.4 Rod Element under Uniformly Distributed Loads 



The uniformly distributed loads are given 
q* ( y . t ) = qo h ( t ) 

Substituting the loads in equation (2.26), we may get the 
consistent nodal forces. 







qo h ( t ) 


l' 


Z 

II 




’ = [Lo + u - u, ] < 








2 


/ 1 " 



du 

where, h(t) = 

dx 

Another case is the linearly distributed 





loads . 



(2.27) 



Figure 2.5 Rod Element under Linearly Distributed Loads 
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Linearly distributed loads are given 



q ( X , t ) = q (t) + [q (t) - q (t)] 




Substituting the loads in equation (2.26) we may get the 
consistent nodal forces. 



D. SOLUTION METHODS 

In this section we derive the equilibrium equations of 
incremental method, Full Newton-Raphson , and Modified 
Newton-Raphson Method. 

It is based on total Lagrangian Method and use the 
Taylor Serier at a given equilibrium point. 

The element stiffness matrices and consistent load 
vectors are assembled to form global matrices in the usual 
manner. The prescribed displacements are implemented by 

replacing the diagonal element of the corresponding row by 
unity and the rest of the element of the rows and columns by 
zeros [ 49 ] . 

1 . Incremental Method (INC) 




(2.28) 



Consider the incremental process, we can rewrite the 



equation of load-displacement relationship, 



[K]i (A>i + i = { A P}i + i + { A Pc }i 



(2.29) 



where, i : the current configuration coordinates. 

{ A Pc }i : the force unbalance 
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{AP}i +1 : an increment in the external load 
The method is described by the Figure 2.6 




where, ABC : The algorithm of equation (2.29) 

AB’ C’ : The algorithm of purely incremental 
process 

Figure 2.6. The Algorithm of Incremental process 

The approximate values of A. E, C, is computed from p =r(u). 
To improve these values, take the iteration process at each 
load level with {AP}i+i = 0, as is Figure 2.7. 

The iteration processes are two ways, one is few 
load levels with many corrective iterations in each, or 
another is many load levels with few iterations in each. 
This is sometimes called incremental with one-step Newton- 
Raphson method". 

2 . Full Newton-Raohson Method (FNR) 

Consider the load versus displacement function P = 
f(u). Take the Taylor series of P = f(u) 
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A,) - + (^-) A A, 



(2.30) 



dP 

where, - = K = tangent stiffness at A 

du A 

A 



f ( u ) = P 

A A 

f (U + A ) = P 

A I B 



We look for A. such that f(u + ^ ) 

’ A ' 



P , so 
B 



K (A) = P - P 

A \ B B 



(2.31) 



The consequent steps are explained by the Figure 2.7. 




Figure 2.7. The Newton Raphson solution of the equation 
P = f (u ) 



The above algorithms are below; 

update displacements : ui = Ua + A , 
use ui to get stiffness ki and load Pi 
. find next increment A s from ki ( A x ) = Pb - Pi 

Eventually ua + A ( + A a + = ub to a close 

approximation . 

Then, the load can be increased from Pb to Pc, and the 

process can be started again. This solution algorithm is 
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known as the Newton-Raphson Method. 



3 . Modified Newton-Raphson Method (MNR) 

The Newton-Raphson method requires that tangent 
stiffness matrix [k] be generated and reduced for equation 
solving in every iteration cycle. This process is not 
economic. An alternative is modified Newton-Raphson 
iteration, also called "constant stiffness iteration" . Here 
the same tangent stiffness matrix is used for several 
iterations cycles. The MNR is described by the Figure 2.8 




As compared with FNR, the MNR requires more iterations, 
but each is done more quickly because [K] is formed and 
reduced in only the first iteration. 

E. ANALYTIC SOLUTIONS 

At first we define the concept of stresses [52] Eulerian 
stress, a , is defined as follows: 
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N(y, t) 

a (y,t) = (2.32) 

A (y) 

and is called True stress. Lagrangian stress, T, is defined 
as follows: 

N(y,t) 

T ( x , t ) = (2.33) 

Ao (x) 

and is called Engineering stress. 

Kirchloff stress, S, is defined as follows: 

N(y, t) 

S(x,t) = — (2.34) 

A(x,t) Ao ( x ) 

The relationship between strain r , and stretch is as 

follows : 

r = S ( \* - 1) (2.35) 

du 

A = 1 + (2.36) 

dx 

where, N : Applied force 
A : deformed area 

For the special case of small strain, both u, x and change 
in cross-sectional area are small. In this case, we can 
write 

r » u _ = € 

i x 

a * T **> S 

where, 6 is engineering strain. 

Thus in small deformation case, the Lagrangian stress T 
and engineering strain G can be used to replace the 
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and engineering strain € can be used to replace the 
Kirchhoff stress S and Lagrangian strain r . 

1 . Elastic Materials 

Consider a one-dimensional rod with varying 
cross- sectional area, Ao = Ao(x), in equilibrium under 
distributed external loading. 




(a) (b) 

where, f : external distributed load per unit length 

T : Lagrangian stress within rod 
Figure 2.9. One-Dimensional Rod under Axial Loading 



For static problems and elastic materials, both f and T are 
independent of time, 

i r £ _ -C f , , \ v- J rp _ 

x v W / i. — +. V / 

From figure 2.9(b), the static equilibrium of forces in the 
x-direction is 

( 2 . 31 ) 

E F(x) = -TAo + f Ax ♦ (T + (dT/dx) ax ) (Ao + (dAo / dx) ax ) = 0 

Expanding in Taylor series, 
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(2.38) 



E'or a continuous rod, 



in the limit, 



i + ^TAJ) =o 



(2.39) 



Substitution of equations (2.34), (2.35), and (2.36) in the 



Note that equation (2.40) gives the values of distributed 
force f (per unit length) required to produce the 
displacement field u(x) of a one-dimensional rod once the 
boundary conditions, the stress-strain law, and varying area 
function Ao ( x ) have been specified. Consider the rod with 
constant area Ao , equation (2.40) becomes 



above 



(2.40) 



where, u(x) : displacement w.r.t. x 



u, x : first derivative of u(x) 



u , x x : second derivative of u(x) 



f(x) = -Ao u.xx [(1 + u,x ) **2 ( dS« /dr ) + S«] 



(2.41 ) 



Some special cases of equation (2.41) is considered. One 



case is the specified displacements. 
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where, u(x) = ao + a, x + a a x 2 
u, k = a, + 2 a x x 
u, xi = 2 a » 

Figure 2.10 Test Model of "Inverse Solution" 

Apply the boundary contition u,x(x) =0 at x = Lo 

a 

u(x) = a, x - — — - x 2 (2.42) 

2 Le- 
if a, is specified, we may get the exact solution u(x) and 
the corresponding required distributed load f(x). 
a. Linear displacements field: 
u ( x ) = a , x 

u,i = a , 

U, X X — 0 



From equation (2.42), we may get 
f (x) = 0 

Note that a linear displacement field cannot be produced 
with distributed load. 
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b. Quadratic displacement field - . 
u(x) = a,x - (ai/2Lo)x 2 
u, x = a, - y ( a « / Lo ) x 

it, xx = - ( a | / Lo) 

- Consider the linear materials, S e = E r 

We may get the distributed forces from equation (2.41) 
f(x) = - EAou, xx [(1 + u. x ) + (u, x + ^ u 2 , x ) ] (2.43) 

- - Consider the nonlinear materials, 

S e = El r + E2 r 2 (quadratic form) 

From equation (2.41), we may get 

f ( x ) = - Ao u, x x { (1 +-u,x) [Ei + 2 E2 (u,x + 'S u 2 , x ) ] 

+ El ( u, x *• *£u 2 , x ) + E 2 ( u, x + ^u 2 , x ) ] } (2.44) 

Another case is the prescribed loads. 




Figure 2.11 Test Model of Rod with Prescribed End Element. 

From equation (2.39), since no distributed loads exist, 
d 

(TAo) - 0 with constant area. 

dx 

We may get 



S 



constant 



(2.45) 



( 2 . 46 ) 



The applied Load is given 
P = Ac \ G 

- Linear material; 

S e - Er = E(u, x + u 2 , k ) 

From equation (2.46), we may get 

P = EAo (u, X + u2,x + u3,x) (2.47) 

- Nonlinear material; 

S e = Ei r + E2 r 2 (quadratic form) 

From equation (2.46) we may get 

P = Ao (l+u,x) [El (u, * +^u 2 ,x)+E2 (u, i + ^u 2 ,x) 2 ] 

( 2 . 48 ) 

2. Quasilinear Vi sc oelas tic materials 

The Viscoelastic materials has the more complicated 
problems than the elastic materials, since the external 
loads and displacements are time-dependent. 



<e — 



U.c-t) 






Ul„C*> 
* 



X=Lo 
— >1 



where , Ao , Lo 
uo ( t ) 
UL o ( t ) 
Fo (t) 
Fl o ( t ) 



Underformed Area and Length 
time-varying left-end displacement 
time-varying right-end displacement 
time-varying left-end force 
time-varying right-end force 



Figure 2.12 One-Dimensional Rod Subjected to Time-varying 
End Loads 
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For quasilinear viscoelastic materials, 



S(t)=S'(t) 



(‘SGU-Q 

I o * 



S e (Qdf, 



(2.15) 



The Lagrangian stress is constant throughout the rod since 
the distributed loads is not applied. From equation (2.39), 
the stress is time variable only. 

T ( x , t ) - g(t) = \(t) S(t) (2.49) 

And the internal forces produced with the rod is 

F(t) = TAo = Aog(t) (2.50) 

Substituting equation (2.15) in equation (2.49) yields 
g(t) = \(t) S(t) 



= (1 + U -x ) 



S e (r) 



‘ dG(t-i) 
o 



S‘(r)dl, 



(2.51 ) 



The right hand side of equation (2.51) is a time function 
only . 

Therefore g(t) must also be a time function, 
du 

= u,x(x,t)=h(t) (2.52) 

dx 

Integrating w.r.t. x at equation (2.52) 

u = u ( x , t ) = xh ( t ) + C ( t ) (2.53) 

where, C : integrating function to be determined by the 

boundary condition at one end. 

Assume that 

u ( o , t ) = uo ( t ) = C ( t ) , 
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therefore 



u(x,t)=xh(tj+uo(t), (2.54) 

u , x ( x , t ) = h(t) 

Substituting equation (2.54) in equation (2.51) yields 

l | t (vfcfelH + 3. U ~S§~ 55) 

Once the G(t) and S e (r) have been specified, the equation 
(2.55) can be used with the boundary conditions to find the 
exact solution in some cases. Consider the 3 element solid 
with linear material property. The prescribed displacements 
boundary conditions at both ends is described in Figure 
2.13. 




B . C . ’ S •• u ( o , t ) = uo ( t ) - 0 

u(Lo,t)=u,(t) (2.56) 

Figure 2.13 Test Model of a Rod with Prescribed Time Varying 
Displacement at the End. 

The prescribed displacements are linear, 
u ( x , t ) - ( x/Lo ) u, ( t ) 

L o 

u, x(x,t) = u. (t) / Lo = h(t) 

l ' o 
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(2.57) 

(2.58) 



where, u^tt) = Loh(t) 

Assume linear material property, S e = Er 

And for a 3 element rod, G(t)is defined below 

G ( t ) = (1 - a ) + e-Pt (2.59) 

where, a , f3 ; constants for material properties. 

Note that G(o) = 1 at t = 0 , G(t) = 1 - a at (3=0 

No relaxation occurs and the material exhibits elastic 
response . 

Specifying the prescribed displacement function u , (t). 

O 

a. Ramp Function 

At x = Lo, B. C.’S : u(o,t) = 0 

u(Lo,t)= T u (t) = at (a: constant) 

From equation (2.57) and 2.58), we may get 

u(x,t) = (a/Lo) x t (2.60) 

u, x ( x , t ) = ( a/Lo ) t 

Substituting equations (2.57) - (2.59) in equation (2.55) 

We may get the reaction force, F. 

EAo at a 2 1 2 a a a — (3 1 

F = [1 + ] [{at+ } { 1- ( )}e 

Lo Lo 2Lo (3 (3Lo 

a a a 

{ ((3t - 1) + ( (3 2 1 2 - 2(3 1 + 2)}] (2.61) 

(3 2 (3Lo 

Equation (2.61) describes the reaction force produced at the 

fixed end for the prescribed displacement 

u , ( t ) = at 
L o 
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b. Step Function 



(Relaxation Test) 



At x = Lo , B.C.’S : u ( o , t ) = 0 
u ( Lo , t ) = ( t ) 



where, l(t) : unit step function 
From equation (2.57) and (2.58), we may get 
UL o 

u(x,t) = x l(t) (2.62) 

Lo 

UL o 

u, * = 1 (t) (2.63) 



Lo 

Substituting equations (2.57) - (2.59) and (2.62) in 
equation (2.55), we may get the reaction force, F. 



F(t) 



EAo 

( ) 



UL o 

UL o ( 1 + 



UL o 

) G(t) [1 + ( ) 



Lo 



2Lo 



Lo 



Consider the Kirchhoff stresses. 



From equation (2.34) 



UL o UL o 

S = E ( ) ( 1 + ) G(t) = E r. G ( t ) 

Lo 2Lo 



1 (t)] 
(2.64) 



(2.65) 



UL o UL o 

where, r, = ( 1 + ) : Lagrangian strain due to 

* Lo 2Lo 

a prescribed step displacement function u, l(t) at the rod 
end . 
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III. PROGRAM IMPLEMENTATION 



A. OVERVIEW 

We are interested in static analysis. For static 
analysis, load-displacement relationship is as following: 

{P} = [k] {u} 

where {P} : external loads vector 

[k] : constrained stiffness matrix 

{u> : nodal displacements vector 
If we know {P} or {u>, then we can analyze the static 
problems by using elementary elastic theory. In general, 
the static analysis steps are shown as follows: 



Step 1 


: Enter geometry of material properties, loads, 
et al . 


Step 2 


: Calculate element stiffness matrices. 


Step 3 


: Determine appropriate part of [k] and {P}. 


Step 4 


: If [k] or {P} is appropriate, then repeat step 2 

through step 3 . 


Step 5 


: If [k] or {P} is not appropriate, then apply 

constraints . 


Step 6 


: Determine {u} . 


Step 7 


: Regenerate elements. 


Step 8 


: Calculate element stresses. 
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B. SOLUTION PROCEDURES 



In nonlinear structural problems, to find the solutions 
we use the INC and FNR and MNR. In fact, it is hard to know 
in advance what increments of loads or iteration should be 
used to obtain a good approximation to the solution. The 
interactive method was used in determining the strategy for 
solution in program application. The strategy subroutines 
were developed for elastic or visco-elastic time-dependent 
structural problems. Here the initial state for incremental 
loading process was selected as a null displacement vector. 
Next we consider finite element method for elastic (FEMEL) 
and viscoelastic (FEMEV) programs. IN both programs, 

equilibrium was checked at each stage for solutions obtained 
from other INC or FNR and MNR. 

1 . Algorithm Implementation 

Consider the algorithm of Newton-Raphson Method 
based on total Lagrangian Method. Full Newton-Raphson 

Algorithm is as follows: 

Step 1. Establish local coordinates X-Y for the element at 
hand. 

Step 2. Generate the element stiffness matrix [k] in local 
coordinates so that it operates on the local D.O.F. 
Step 3. Transform [k] to global coordinates so that it 
operates on global D.O.F. 

Step 4. Repeat step 1 through 3 until all elements have been 
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treated. 



Step 5. Compute displacements. 

Step 6. Compute internal forces, [k]{u}. 

Step 7. Solve [k] (Au} = (A p > 
where, u N+1 = u N + A U N 

Step 8. Test for convergence. 

If not satisfied, return to step 1. 

Where, convergence = ((A pT AR) / (Ap t Ap))**0.5 < 
0.001 
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2. Flow Chart of the Fi ni te E lement Program 




piowcha^^TthTTirntT^IemerU^^ 
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I v . NUME RIC AL EXAMPLES 

We consider two numerical examples for elastic and 
viscoelastic materials. 

A. ELASTIC MATERIAL 

1 . Prescribed conce n trat ed forces . 

Specify the rod element. 




P, 



Pa - 22,5000 lbs 
A = 0. 25 [ins ] 

L = 10 [in] 

E - 106 [lbs/ins ] 

Consider the stiffness matrix of rod element, [k] . For 
the rod element, 



EA Ki i ki 2 -i 

[K] = 

L |_k2i k2 2 

Consider the boundary conditions, 
ui = 0 
Pi = 0 



(4.1) 
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p2 = 22,5000 at x = L 

From loads - displacements relationship, 



EA 


r 1 


_1 1 






L 


L-i 


iJ 


1 U2 J 


lP2. 



Solve for the displacements, 
Hence, we may get 



U2 



LP2 

U 2 = = 9.0 ( in ) 

EA 



Pres cr ibed di splacem e nts 




X 



Consider the same rod element before, and specified 
the quadratic displacements, u(x), 
u(x) = 0.02 x - 0.001 x 2 

The first derivative is 
du 

- - 0.02 - 0.002 x 

dx 

The second derivative is 

d 2 u 

- 0.002 

dx 2 
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From equation ( 2.41 ) with the constant area, we may 



getthe distributed loads for any x. 



d 2 u du ds e 

f ( X ) = - Ao [(1 + )2 

dx 2 dx dr 

= 3.0 * 10-3 X 2 - 3.06 x + 



+ S e ] 

530 (4.2) 
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FORCE .VS. NO. OF ELEMENTS 




S33H0i 



4 



0 10.0 1 Z .0 M.0 10.0 10.0 zo.o n.o 

NO. OF ELEMENTS 



DISPLACEMENTS .VS. COORDINATES 




4 3 



0.0 1.0 14 a.o 4.0 e.o o.o 7.0 0.0 0.0 10.0 11.0 

COORDINATES 




STRESSES .VS. COORDINATES 




fiu sassaais iiOHHoaiH 
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COORDINATES 



B. VISCOELASTIC MATERIAL 



Consider the rod element below, 

► U^-Jp) 

EkliO, 

"> X. 






U.c^) 



Boundary conditions; u(0,t) = uo (t) = 0 

u ( Lo , t ) = u (t)=u l(t) 

L o t o 

Lagrangian Stress; T(x,t) = g(t) = A(t) S(t) 
Internal Force; F(t) = T Au g(t) 

P 

Kirchhoff stress; S - -r 

A ( x . t ) Ao ( x ) 

For 3 element solids, the relaxation function 
G ( t ) = (1 - a) + a e (3t. 

where, Oi and 3 are constants of the 



(4.3) 
is 

(4.4) 
material 



properties . 

Assume unit step function for prescribed displacements. 

At that time, the boundary condition is 
u ( o , t ) = 0 

u ( Lo , t ) = u r ( t ) = u. 1 ( t ) 
where, l(t): unit step function 

We may get the displacements function with unit step 

f unction . 



x 

u ( x , t ) - ul o 1 ( t ) (4.5) 

Lo 

4 5 



du 



UL o 



l(t) 



(4.6) 



dx Lo 

And the Kirchhoff stress is specified 

S = E r (4.7) 

Substituting equations (4.4) - (4.7) in equation (4.3), 

EAo UL o UL o 

F(t) = Ulo 1 + G ( t ) 1 + l(t) (4.8) 

Lo 2Lo Lo 

Hence, we may get the Kirchhoff stress, S, 

S = r, G ( t ) (4.9) 

L o 

UL o UL o 

where, rLo = 1 + : Lagrangian strain 

Lo 2Lo 

due to the prescribed step displacement function 
u l(t) at the rod end. 

Equation (4.9) describes the relaxation cases. From 
equation (4.8), we may get the internal forces. 
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RELAXATION .VS. TIME 



9 




o 

b 
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TIME 



V. CONCLUSION 



This study is directed towards understanding the 
nonlinear behavior of structures. 

The formulation is based on the principle of virtual 
work. We consider the one dimensional nonlinear continuum- 
based rod element that includes large displacements and 
viscoelastic effects. To get the solution we use the 
incremental process with full Newton-Raphson and modified 
Newton-Raphson Method. Numerical solutions agree well with 
exact solutions. 

Further studies may be extended to two dimensional 
structures, such as laminated composite plates. 
Consideration of the viscoelastic effects of composite 
panels at high temperatures is another possible area of 
study . 

An experimental study to verify the results of the 
present study is also recommended. 
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